Exact mesoscopic correlation functions of the pairing model 
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We study the static correlation functions of the Richardson pairing model (also known as the 
reduced or discrete-state BCS model) in the canonical ensemble. Making use of the Algebraic Bethe 
Ansatz formalism, we obtain exact expressions which are easily evaluated numerically for any value 
of the pairing strength up to large numbers of particles. We provide explicit results at half-filling 
and extensively discuss their finite-size scaling behavior. 



I. INTRODUCTION 

The pairing phenomenon is ubiquitous in quantum 
many-body systems of sizes ranging from the very small, 
like quarks and nuclei, to the very large, like stars 1 ^. The 
common feature of all these seemingly unrelated systems 
is the instability against the formation of Cooper pairs for 
an arbitrarily weak attractive force, the basis of the BCS 
theory of superconductivity^. Despite the diverse nature 
of pairing systems, many of their fundamental properties 
can be understood phcnomenologically from a so-called 
reduced BCS model. 
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which was introduced by Richardson in the early 1960's 
in the context of nuclear physics^. The model simply de- 
scribes (pseudo) spin-1/2 fermions (electrons, nucleons, 
etc. . . ) in a shell of doubly degenerate single particle 
energy levels with energies e a /2, a — 1, ... N. c a ,a are 
the annihilation operators, a = +, — labels the degen- 
erate time reversed states (i.e. spin or isospin) and g 
denotes the pairing coupling constant. Despite its simpli- 
fied character (all levels interact uniformly), this Hamil- 
tonian captures the main essence of the problem; fixing 
(ab-initio or phcnomenologically) the energy levels e a and 
the coupling g allows to obtain quantitative predictions, 
since the model remains solvable for an arbitrary choice 
of parameters. 

In the thermodynamic limit, and within the grand- 
canonical ensemble, the properties of the Richardson 
model are correctly described via the BCS variational 
ansatz 3 . However, for finite numbers of particles, the 
situation is more complex. The actual solution then de- 
pends on the ensemble chosen, and for physically rele- 
vant systems the grand-canonical ensemble is not always 
the appropriate one. For example, nuclei have a fixed 
number of nucleons; due to the typically large charg- 
ing energy, experiments on ultra-small metallic grains 
are also performed at fixed number of electrons 5 . In 
these cases a treatment based on the canonical ensemble 
would be more appropriate, precluding a BCS mean-field 
approach. Moreover, dealing with a system in the meso- 
scopic regime precludes the use of quantum statistical 



mechanics, and one is thus forced to rely either on un- 
controlled approximations, or nonperturbative methods. 

Fortunately, Richardson's Hamiltonian ([I]) is one of the 
theories for which an exact solution can be constructed in 
the canonical ensemble 4 . This solution explains several 
interesting features of the mesoscopic physics of super- 
conductors, complementing previous approximate treat- 
ments (see the review jf|). In particular it allowed to 
give a definitive answer to Anderson's 1959 question^: 
What is the size limit for a metallic grain to have su- 
perconducting properties? The utility of the model is 
thus indisputable (see also the reviews 8,9] for some 
non condensed-matter applications), however most of 
the attention has been concentrated on thermodynamical 
quantities. On the other hand, experiments typically give 
access to static or dynamical correlation functions, which 
are not easily obtained in this framework. Richardson 
himself in 1965— derived a first exact expression for static 
correlation functions, which unfortunately has a degree of 
complexity that grows factorially with system size, and 
was therefore not suitable for actual calculations. In a 
significant development, Amico and Osterloh 11 proposed 
a new method (based on a generalization of earlier work 
by Sklyanin 1 ^) to write down such correlations explicitly. 
The complexity of this method was still factorial and all 
the numerical results were therefore limited to system 
sizes of up to 16 particles. A disadvantage of these meth- 
ods is that all the eigenstates of the Hamiltonian must 
be known to get the correlation functions. 

A major simplification was then proposed by Zhou et 
ali 1 ^ (see also [lj])- Using the Algebraic Bethe Ansatz 
(ABA) and the Slavnov formula for scalar products of 
states^ 5 -, they managed to write the static correlation 
functions as sums over N$ determinants of N r x N r matri- 
ces, reducing the complexity of the problem to order 
(here N r is the number of rapidities in the eigenstates). 
Furthermore, in this approach only the knowledge of the 
ground-state wavefunction is required. Surprisingly, this 
approach has not been used until now to obtain quantita- 
tive numerical results for the correlation functions, with 
the notable exception of the calculation of ground-state 
entanglement properties 1 ^. 

In this paper we fill this g ap. As a first step we re- 
analyze the results of Refs. [13lll4| . rewriting all static 
correlation functions as sums over only N r determinants 
(thereby reducing the complexity of the problem by a 
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further factor of N r ). We then provide analytical formu- 
las for the physically relevant correlation functions, and 
evaluate them for some model Hamiltonians. We stress 
that having reduced the complexity of the problem by 
this amount, correlation functions of systems with many 
more particles than before can be calculated on a sim- 
ple computer. In this way we can describe the crossover 
from mesoscopic to macroscopic physics, going beyond 
previous results limited to fewer particle o 11 ' 17 . 

The paper is organized as follows. In Sec. [II] we dis- 
cuss the model and its general properties. In Sec. Unl 
we recall how to calculate static correlation functions by 
means of Algebraic Bethe Ansatz, and derive their gen- 
eral expressions in terms of sums of N r determinants. 
This section is rather technical, thus the reader inter- 
ested in the physical results can skip directly to Sec. IIVI 
where we discuss how to solve the Richardson equations 
for the ground state, and derive quantities that do not 
require knowledge of the determinant representation. In 
Sec. |V]all the correlation functions are explicitely calcu- 
lated at half-filling. The paper is closed by Sec. I VII where 
we also discuss open problems for future investigation. 



II. THE MODEL 

A simple but very important property of the system 
is the so-called blocking effeclj 4 ^, i.e. unpaired particles 
completely decouple from the dynamics and behave as 
if they were free. We will denote the total number of 
fermions as Nf, and the total number of pairs as N p . 
Due to level blocking, we will only consider Nf — 2N p 
paired particles in N unblocked levels. In terms of pair 
annihilation and creation operators 
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the Hamiltonian is 
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and n a — 2b' i a b a is the number of particles in level a. 

The pair creation and annihilation operators satisfy 
the commutation relations 

K, ftj] = <Mi - , [b a , bp] = [ftt , ftt] = o . 

(4) 

The term 2b\,b a in the first commutator makes the model 
different from free bosons and therefore non-trivial. 

Using the pseudo-spin realization of electron pairs 
S z a = blb a - 1/2, S~ = b a , S+ = bl the BCS Hamilto- 
nian becomes (up to a constant) 
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The operators S^' z obey a standard spin algebra and so 
the Hamiltonian ([5]) describes a spin- 1/2 magnet with 



long-range interaction for the XY components in a site- 
dependent transverse magnetic field e Q . Such a mag- 
netic Hamiltonian is known in the literature as a Gaudin 
magnet 18 . An important relation is 



Q± QT — Q 2 _ I Q z \ 2 -L- Q z 
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Since the normalization of the pairing strength in the 
literature is not uniform, care must be taken when com- 
paring the results we will obtain with other papers (e.g. 
our g is the half of the coupling used in Refs. HIII3])- 



A. Grand-canonical BCS wavefunction 

In the grand-canonical (GC) ensemble the ground- 
state wavefunction is the BCS variational ansatz 



\GS) =H(u a +e l ^v a bi)\0), with^ 



vl = 1. 



(7) 

where the variational parameters u a and v a are real and 
(f> a is a phase which, it turns out, must be a- independent. 
\GS) is not an eigenstate of the particle number operator 
Nf and the average condition (Nf) — Nf determines the 
GC chemical potential. Likewise, the commonly used 
definition 



A GC = 2 g y(b a ) 
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for the superconducting gap makes sense only in a GC 
ensemble, since (b a ) is zero when evaluated at fixed par- 
ticle number. The variational parameters are obtained 
as 



y/(e a - A*) 2 + |A GC f 



(9) 



where /i is the GC chemical potential. 

It is then easy to calculate (static) correlation functions 
on this GS: 
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B. Canonical description and Richardson solution 

The exact solution of |T]) was worked out by 
Richardson-. In the canonical ensemble the model 
is integrable^ and tractable by means of algebraic 
method o 13 i 14 ' 20 i 21 . We review here only the main points 
of this solution. 

In the ABA, eigenstates are contructed by applying 
raising operators on a so-called reference state (pseu- 
dovacuum). We here choose the pseudovacuum (in the 
spin representation) to be fully polarized along the z axis 
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In the pair representation, this state thus corresponds to 
having one pair in each available level. Eigenstates with 
N p pairs are then characterized by N r = N — N p spectral 
parameters (rapidities) Wj, and take the form of Bethe 
wavefunctions 
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The operators B, together with operators A, C, V defined 
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obey the Gaudin algebra, which is the quasi-classical 
limit of the quadratic Yang-Baxter algebra associated to 
the gl(2) invariant i?-matrix (we refer the readers to [3| 
for details). 

The wavefunctions (fT2|) are eigenstates of the trans- 
fer matrix, and thus of the Hamiltonian |T]), when the 
parameters Wj satisfy the Richardson equations 
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Throughout the paper we will refer with latin indices to 
the rapidities and with greek ones to the energy levels. 
The total energy of a Bethe state is E = J2a=i 1? — 



- g(2N r — N). For a given N and N r the number 



of solutions of Richardson equations is ( N ) , and coin- 
cides with the dimension of the Hilbert space of N r pair 
vacancies distributed into N different levels, i.e. the so- 
lutions to Richardson equations give all the eigenstates 
of the model. 

Note that the Richardson equations fT4")) have a differ- 
ent sign of g compared to the ones mostly considered in 
the literature. This is due to the particular choice of the 
pseudovacuum we made in Eq. (TTj) , whereas the most 
common choice is S^\0) = — 1/2 10). With our choice 
of pseudovacuum, Bethe states are built by destroying 
pairs, as in Eq. (fT2")l and not by creating them. We 
use this somehow unusual pseudovacuum following Refs. 
[l3lfl4j in order to use all the formulas there without any 
adaptation. At half-filling (that is the only case consid- 
ered numerically here) the different choice of the pseu- 
dovacuum only matters as a global normalization and a 
different labeling of the states. 

The connection between the canonical and grand- 
canonical ensembles was first pointed out by Richardson 
himself 22 , who showed how in the large Nf limit one re- 
covers the BCS gap equation as 
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where now fi is fixed by the density and the equation can 
be solved to find A, which with this normalization is an 
intensive quantity and corresponds to Aqc /Nf . For the 
ground-state energy per pair Eq one finds 
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Anderson^ argued that increasing the mean energy spac- 
ing d (that is inversely proportional to the volume in a 
metallic grain) superconductivity should disappear when 
d becomes of the order of the bulk gap Aqc ■ Our study of 
correlation functions to be presented below clearly shows 
this crossover. 



III. ALGEBRAIC BETHE ANSATZ AND 
CORRELATION FUNCTIONS 

The starting point to calculate correlation functions 
with the Algebraic Bethe Ansatz (ABA) is having a rep- 
resentation for the scalar products of two generic states 
defined by N r rapidities (N — N r pairs) 
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({w}\{v}) = {0\l[C(w b )l[B(v a )\0), (17) 
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when at least one set of parameters (e.g. w b but not 
v a ) is a solution to the Richardson equations. Following 
standard notations, C is the conjugate of the operator B. 
Such a representation exists, and is known as the Slavnov 
formulaic, which for the case at hand specifically readsi^ 
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Ub<a( W b ~ W a ) Ua<b( V b ~ Va) 

xdet Nr J({v a },{w b }), (18) 



where the matrix elements of J are given by 
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from which the norms of states simply follow from v — > w 
as Hl^ill 2 = det/v.. G with a Gaudin matrix 



N Nr 

E ? ^ - 2 E?-^ 

( ■»)_ — fo V ill- — 



c^a 



(V a - V c ) 2 



a = b , 

a ^ b, 
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recovering Richardson's old result^. 

The key point is that any form factor of a local spin op- 
erator between two Bethe eigenstates can be represented 
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via (|T3)) as a scalar product with one set, e.g. {v} not 
satisfying the Bethe equations, for which Slavnov's for- 
mula is app licable. This has been explicitly worked out 
in Ref. [lj]. For {«;},{?;} containing respectively N r + 1 
and N r elements, the nonzero form factors arc: 



nlV r +l/ 
6=1 [Wb-f-c 
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detN T +iT(a, {w}, {v}) 
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and, for both {w} and {v} containing N r rapidities 
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with the matrix elements of T given by 
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Above, T z is the N r x AT r matrix obtained from T by 
deleting the last row and column and replacing N r + 1 
by N r in the matrix elements. Here it is assumed that 
both {v a } and {wb} are solutions to Richardson's Bethe 
equations. However, the results are still valid for S± if 
only {wb} satisfy the Bethe equations. 



A. Determinant representation of the correlation 
functions 



In Ref. [13| it has been pointed out that due to the 
simplicity of the solution of the ABA not only the form 
factors, but any static correlation function can be written 
in a determinant representation. This simplicity puts 
the BCS model in an extremely privileged position for a 
detailed study of the static correlation functions. 

The result for ({w}\S~ Sp\{v}) has been explicitly 
worked outi^ 

({w}\S-S+\{v}) = ]T ^—{{ W }\s-\{v}i) 
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Here the sets indicated by {v}i stands for sets where the 
rapidity i has been removed and similarly for both 



i and i' rapidities have been removed. The S S form 
factor is given byi^ 
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where a ^ (3 is assumed, with the convention that it 
vanishes when a = j3. Note that ({w}|S'~S'J|{w}) is 
symmetric under the exchange of a and /3, although this 
is not manifest in the formal expression. This nontrivial 
property will be checked during the numerical computa- 
tion. 

This correlation is then written as the sum of N% de- 
terminants, which is much less than the sum over the 
full Hilbert space needed in other approaches. In the fol- 
lowing we will determine a similar expression for (S^S^) 
and then we will show that it is possible to reduce these 
formulas to sums of only N r terms. 



1. Determinant representation of (S^Sp) 

The operator A{u) only has simple poles at e a such 
thaOi 



S* — lim (u — e a )A(u) 

u — *e a 

This allows one to write 



{{w}\S* a S* p \{v})= lim lim ( u ' - e a )(u - e p ) 
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which can be easily computed by commuting the A oper- 
ators until they reach the far right and act on the pseu- 
dovacuum in the following way: 
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Using the commutation relation^ 

B(u) B(v) 
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and denning Q = ({w}\ A(u')A(u) Jl -Ji |0), we find 
by commuting A(u) and B(v\) that 
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where |{f }i; u) is the (non-Bethe) state built by replacing 
the rapidity v\ by u. By commuting again N r — 1 times 
and using Eq. (j2"9"l) , we find 
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where we defined F(u) 



The same procedure can then be repeated in order to 
have A(u') act on |0): 
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with Fi(u) e-I + IV", -jl -r^. It is 

then easy to take the limit as prescribed by Eq. (|28|) and 
find 
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Note that we use the same notation as in Eq. (f2Tjl for 
the S~ form factor but there it is evaluated between two 
states with N r and N r + 1 rapidities, while here it is 
between two states with N r — 1 and N r . This should not 
be a source of confusion. 

The static correlation function can now be evaluated 
by setting {v} = {w} = {wo}; the set of rapidities corre- 
sponding to the ground state of the system. Using equa- 
tions (|2"Tj) and (|2"5")) , we can directly express this correla- 
tion function as a sum of iVj? + N r determinants. 



B. Reduction formulas 

As anticipated in the introduction, the previous ex- 
pressions can be reduced to sums over only N r determi- 
nants. This is explicitely worked out in the following. 
We will assume that a ^ (3, because for intra-level cor- 
relations from Eq. ^ we have 

(Ml S-Si \{w}) = \ + (M| SI \{w}) , (35) 
which is already a single determinant expression. 



1. Reduction of (S a St) 



We here need to evaluate Eq. ([24]) in the limit v — > w. 
In this case, Eqs. (|2"Tj) and (|25|) for the form factors 
simplify to 

(MaU^IWaw) = ™«-*° det „ (36) 

w q - ep w q -ep 
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The matrices U are defined as follows. and f/(''Q 

are equal to the Gaudin matrix (|20p except for columns 
q and q, I respectively, where 
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In Eq. (|37|) it is explicitely assumed that a ^ (3 with 
the convention that for a = (3 it is zero. And in fact, for 
a = /3 it is not difficult to prove that Eq. (f31)]) reproduces 
the correct result given by Eq. (|3"B")) . 

Since everything is symmetric under exchange of I and 
g we only perform the sum over I < q and in the end 
we multiply the result by 2. Thus we need to perform 
the sum (neglecting for the moment the / independent 
factors) 



y ^LZI^. dct = y K lq det 

f-^Wl-Wq ^ ' 



1 = 1 



(40) 



Let us write the matrix U^ ql > as a vector of vectors 

jjiw = |g 1 _ _ _ Gi-i,B, G/+i . . . Gq-i,C, Gq+i ■ ■ ■ GjvJ , 

(41) 

where the Gi corresponds to the columns of the matrix 
G, C (at position q) corresponds to the vector given by 
Eq. O and B by Eq. 
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The sum we want to calculate is (we use 
determinant) 



for the 
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K^detuW = K lq \B,G 2 ,G 3 .. 



1=1 



K 2q \G 1 ,B,G 3 ...\+K 3q \G 1 ,G 2 ,B 



(42) 



Using the fact that two determinants which differ by 
a single column can easily be expressed as a single de- 
terminant the two first terms of the sum can be writ- 
ten as \K 2q Gi — Ki q G 2 ,B, G 3 . . . |. Elementary column 
operations allow us to write the third determinant as 

<2, 



K^q \Gx, G 2 ,B . 



\G X 



K2 ':G 2 , G 2 ,B . 



This 



l3 ? |Lxl, Lx2, £> . . . I = K 3q \ 

term differs by a single column from the preceding sum. 
We can then write the sum of the first three terms as a 
single determinant K 3q \G x - ^G 2 , G 2 - #^G 3 , B . . . I. 
We can keep on adding terms in the same way until we 
reach column q — 1 and find that Eq. (|42|) can be written 
as a single determinant 

K q -i q \G\ — ——^-G 2 , G 2 — —^-G 3l . . . B,C, Gq+i ■ ■ ■ GN r \ ■ 



2<; 



(43) 

In this way we reduced the double sum to a single one. 
The additional terms in the correlation function (coming 
from (S~)) can also be incorporated to the sum in a 
similar fashion. The (S~) term in Eq. (|24|) is given by 
Eq. ()36|) and can be simply encoded in the representation 
we just obtained for the sum over I < q of (S~S~). In 
this way, we finally have for the full correlation function 

({w}\S-S+\ W) = £ ?±^DM , (44) 



2. Reduction of (SaSp) 
According to Eq. ([34]). when v — ► w 



({w}\S*S%\{w}) = \({w}\{w}) 
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~ Wi - e a (wi - ep)(wi> - e a ) 



(47) 



^From Eqs. (|44l) and (|45| it is straightforward to show 
that when a ^ (3 



({w}\S^\{w}) = 

ll ll 2 1 Nr 

il^-i^(det^)+det^)). (48) 

9=1 

For a — f3 the result is trivially ({w}|(S^) 2 |{u;}) = 1/4. 

This completes the representation of the static correla- 
tion functions in terms of determinants. To make further 
progress, we need explicit results for the ground-state 
rapidities {w}, i.e. the lowest-energy solutions to the 
Richardson equations. The following section is devoted 
to this. 



IV. THE SOLUTION OF THE RICHARDSON 
EQUATIONS FOR THE GROUND STATE 

A. General properties 



where we defined the matrix 
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that has the following structure 
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The low level of complexity of this representation as sum 
of iV r determinants of N r by N r matrices allows us to 
access easily the static correlation functions for systems 
with large number of pairs compared to previously pub- 
lished results. 

We stress again that these formulas are true only for 
a ^[3. 



At g = 0, for pair vacancies in N energy levels e a 
with only double degeneracy, the ( N ) solutions to the 
Richardson equations are trivial. They are given by Eq. 
(fT2|) with the N r rapidities set to be strictly equal to 
one of the energies e Q . Clearly, the GS in that limit is 
built by choosing the N r highest energy levels, i.e. w\ = 
£jv i w 2 — ejv-i ■ ■ ■ WN r — ejv-Ar r +i- 

Apart from a few particular cases with a small number 
of particles, the Richardson equations are not solvable 
analytically when g ^ 0. The perturbative expansions for 
small 23,24 and larg o 25 ! 26 coupling are not predictive for all 
values of the pairing strength and so the most accurate 
results come from the numerical solution. The solutions 
are such that every Wj is either a real quantity or forms, 
with another parameter , a complex conjugate pair 



(CCP), i.e. w*, 



3 ' 



The mechanism for the CCPs 



formation is very easy: as interactions are turned on, 
all Wj are real quantities for small enough <?, but at a 
certain critical value of the coupling g* two rapidities 
will be exactly equal to one of the energy levels (wj = 
Wj' = e 7 (j)) and for g > g* , the two parameters that 
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FIG. 1: g dependence of the real (left) and imaginary parts (right) of the ground state rapidities. From top to bottom they 
correspond to N r = 8, 15, 64 always at half-filling (N = 2N r ). 



collapsed will form a CCP at least for a finite interval in 
g. The situation is in fact rather intricate: the values g* 
are implicit functions of all other rapidities, and can only 
be read off a full solution of the Richardson equations 
for a specific choice of state. Moreover, CCPs can split 
back into real pairs, whose components can then re-pair 
with neighbouring rapidities. Finding complex solutions 
to the Richardson equations is thus difficult in general, 
since there is no equivalent to the 'string hypothesis' as 
for e.g. integrable spin chains. 

The solutions for the ground state have a particularly 
simple structure, fn fact, the set of critical points is such 



that the smaller a rapidity is at g — the smaller the g 
at which it forms a CCP will be. As we raise g from zero 
there will come a point at which u>jv r will form a CCP 
with WN r -i when they are both equal to eN-N r +2- As g 
is raised some more, wm t -2 and w^ r -3 will form a CCP 
at £N-N r +4 and this will go on until every rapidity has 
formed a CCP in the case of even N. Oppositely, with an 
odd number of rapidities in the system, w\ (the largest 
one at g — 0) will always remain a real quantity no matter 
how large the coupling strength is. After the CCPs are 
formed no further collapse happens in the case of the 
ground state, while for excited states further collapses 



8 



can take place and complex solutions can become real 
again. 

Different choices of the parameters e a and of their 
eventual degenerations specify different models. In all 
the preceding sections everything was completely gen- 
eral (modulo having to take some extra precautions in 
the case of coinciding levels e Q ), but from now on we 
specialize to the case of equally spaced doubly degener- 
ate levels. We make the choice to use e a = a which 
sets the zero of energy and implies that every energy will 
be given in units of the (pair) inter-level spacing. Fur- 
thermore we consider only half-filling of the energy levels 
(N = 2N r = 2N P = N f ). In this case, as g — > oo, the 
real part of every rapidity will go to +00 whereas the 
CCPs imaginary parts will go to ±00. 



B. Numerical procedure and results 




70 80 90 100 110 120 130 140 150 160 
Re[wj] 

FIG. 2: Location in the complex plane of the ground-state 
rapidities showing the formation of the arc-like solution. All 
the values correspond to N = 128, N r = 64. 



tions 



At the precise value of g at which a pair of rapidities 
(Wj , Wf ) collapse into a CCP (wj = Wj> — e 7 (j)), the 
Richardson equations (Eq. (|14|)) labelled j and f will 
include two diverging terms whose sum remains finite. 
In order to be able to treat these points numerically, one 
can define the following real variables, 



Wl j = Wj 



Wj' 



2e 7 (j) - w 3 



(Wj-Wj,) 

whose inverse transformation reads 



Wj 



(49) 
(50) 



1 



in-. 



' 2ej_i - wij 

W2,j 



' 2ej_i — wi, 

W 2,j 



(51) 
(52) 



As discussed in Ref. [23|, we need to know beforehand 
which rapidities will form a CCP and at which e 7 (i) it 
will happen in order to use this type of change of vari- 
ables. Since in this article we only need ground state 
solutions, this requirement is easily met. 

At the critical point W2j goes to a well defined (though 
a priori unknown) finite 0/0 form. Using it as a variable 
in the system of equations therefore avoids some potential 
numerical complications when close to a critical point. 

By multiplying the j and j' Richardson equations re- 



spectively by e 7 (j) 



and e 7 (j) — Wj>, we can get 



rid of the divergences that show up at critical points. 
Adding the resulting equations (giving Fij) and sub- 
tracting them and then dividing it by Wj — Wj> (giving 
i<2j), it is simple to obtain the following two real equa- 



^ = E 



l f -12 _ ( F . _ m*L) w , - - ^j-i-WLj 



N, 

E ^ 



(ej_x - w r )(2w f - WlJ ) 



1 2j 



2e i _i - wij 



N 

E 



+ 2(AT-l)-4(AT r -2) = 0, (53) 



2 (e 3 _i - Wf) 



1 

9 



+ 2w 2 ,j = 0. (54) 



The resulting system of non-linear equations can then 
easily be solved using Newton's method. Notice that 
every element of the Jacobian matrix has an analytical 
expression that is easy to obtain and therefore is not ex- 
plicitly written here. Of course, for Newton's procedure 
to converge to the correct solution at a given g, we need 
a good approximation to it. It is simple to do so by 
slowly incrementing g starting from g — 0, where the GS 
is known. One can then use a simple linear regression 
on wxj,W2,j to obtain an educated guess to the ground 
state at g 1 = g + Ag. Despite its simplicity this method, 
very similar to other ones in the literatu re 6 ' 23 ! 27 ' 28 , is 
sufficient for obtaining the ground state solutions. For 
general states, for which the formation of and splitting 
apart of CCPs can be highly non-trivial, a more refined 
algorithm (see Refs. [29ll30l | for example) is needed to 
find the solutions. 



9 



0.6 
0.5 
0.4 
0.3 
0.2 
0.1 


















O 


o 


o 










rf 

! 
! 

jo 












.._ 


■ 






First CCP 
Last CCP 
[2arcsinh(1)]~ 


• 

o 




- 




















- 


• 
• 
• 


















V 
















• 


•- 




• 





50 


100 


150 200 


250 



FIG. 3: Values of the coupling parameter at which the first 
and the last complex conjugate pairs (CCP) are formed 
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FIG. 4: Ground state energy per pair Eq as a function of g. 
Inset: Zoom close to the crossing point 



Fig. [T] shows three examples of the numerically com- 
puted ground state solution of Richardson's equations. 
One can see that the generic statements made about this 
solution in the preceding subsection are confirmed. As 
g gets sufficiently large and every rapidity has collapsed 
into a CCP (for an even number of pairs), they arrange 
themselves into an arc in the complex plane as shown 
more clearly in Fig. O This behavior was originally pre- 
dicted using the analogy between the set of equations and 
a two dimensional electrostatic problem involving fixed 
and free charge o 18 i 22 ' 31 . 

For a correct interpretation of the main features of the 
solutions to the Richardson equations it is important to 
know the value of the superconducting gap given by Eq. 
(fT5)) for the particular Hamiltonian we choose (i.e. e a — 
a). It is easy to show that for large N 



±GC 

N 



1 



2sinhl/2g ' 



(55) 



an expression we will need to compare finite-size results 
with the grand-canonical ones. Consequently Anderson's 
criterion^ for the presence of superconductivity for large 
N is 



A > N 



> 



2 In 2N r 



(56) 



showing the typical 2 ^ logarithmic behavior of the small 
g expansion. 

Fig. [3] shows, as a function of N r , the values of the 
coupling constant g* N (N r ) at which the first two rapidi- 
ties form a CCP. We also plot, for even N r , the values 
of g = g*(N r ) at which the last couple of rapidities col- 
lapses into a CCP. The latter is limited at large N r by^L 
go = (2arcsinhl) _1 = 0.567296 a constant which is also 
shown in the figure. 

These two numbers are particularly relevant to under- 
stand qualitatively the different behaviors as function of 
g and N r . In fact, when g is larger than g*(N r ) all the 



particles are paired and the system has entered its asymp- 
totic superconducting regime. Oppositely when no pair 
has still been formed, i.e. for g < g% (N r ) supercon- 
ductivity is absent. In fact, g* N (N r ) coincides with the 
critical value of the coupling given by Anderson's crite- 
rion Eq. (|56[) for large N. The curve resulting from Eq. 
(1S"6|) is plotted in Fig. [3] and the agreement with g^f r {N r ) 
is excellent even for relatively small value of N r . 

Note also that g* N (N r ) vanishes in the thermodynamic 
limit, which can simply be interpreted as the Cooper 
instability. A quantitative understanding of these phe- 
nomena and of the crossover between small and large g 
at fixed finite N r requires an accurate study of the cor- 
relation functions, which we present in the next section. 



C. Ground state energy 

In Fig. [4] we plot the value of the ground state energy 
per pair (in units of the inter-level spacing) at half-filling 
for a set of different number of pairs as given by N p Eq = 

One interesting feature is the presence of a size invari- 
ant point at which every curve cross (see the inset of Fig. 
|4]for a zoom close to this point). Indeed at gi nv w 0.910 
the ground state energy seems to be independent of the 
number of pairs in the system E(g inv ) « —0.45. How- 
ever, the presence of this "fixed point" does not carry any 
deep meaning and can be easily understood in terms of 
the 1/Np expansion developed in Refs. [22j|32|. In fact, 
according to these references for large N p the ground- 
state energy per particle can be written as 



En 



N p E 



(o) 



E 



(i) 



0(1/N P ) 



(57) 
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FIG. 5: (Color online) function of a. Each plot is at 

fixed g and for several N p = N/2. 



with [E^ 0) is nothing but Eq. JTBJ)] 



E 



(0) 


(1) 



E 



<t>(2g) 



1 coth — , 

2 2 fl ' 

1 



l-0(2.g) coth l/(2.g)), 

dx cosh7ra;/2 



(58) 
(59) 



1 + X' 



^cosh 2 (7rx/2) + sinh 2 (l/(2 5 )) 



where we adapt the results to our normalization (i.e. 
the quantities of Ref. [32[ reads D = 2N p , A = 2.9 
and there is a global shift of the energy levels). The 
scale invariant point just corresponds to the value of 9 
for which the order N p term E^ vanishes, i.e. gi nv = 
(2arccoth2) _1 = 0.910239 .... The energy at this point, 
apart from 0(1/N p ) corrections, is independent of N p 

and given by E^ ] (g mv ) = -0.45276.... Eq. (57|) is 
thus practically a perfect approximation of the actual 
value of the ground-state energy for large enough N p , 
say N p > 16. 



V. CALCULATION OF THE CORRELATION 
FUNCTIONS 

The formulas we obtained for the correlation functions 
are completely general and are valid for any choice of 
the Hamiltonian parameters e Q and g (some care would 
have to be taken in the limit of coinciding energy lev- 
els, however). To obtain a physical result we still have 
to perform the sum over the N r terms, introducing in 
the determinants for the form factors the solution to 
the Richardson equations. This cannot be done ana- 
lytically, so we need to make a choice of the model to 
study. As we already mentionned, we only consider the 



FIG. 6: (Color online) function of a. Each plot is at 

fixed N p but for several different g going from 0.1 (always the 
smallest) to 1 (always the largest) increasing by steps of 0.1. 



most-studied case in the condensed matter literature, 
which consists of N equidistant levels at half-filling, i.e. 



N = 2N r — 2N n 



Nf. We normalize the levels as 

a with a = l...N, (60) 



i.e. we measure the energy scale in terms of the inter- 
level spacing and we fix the Debye frequency (the largest 
energy level) to N. 



A. Correlations among the same level and 
"canonical" order parameter 

Among the various correlation functions a central role 
is played by the ones on the same level. We consider the 
correlation 



u a v D 



(61) 



that can be easily obtained by the previous representa- 
tion of (S^) and does not require the reduction formulas 
because it is written in terms of a single form factor. This 
correlation is important because it is one of the building 
blocks of the BCS theory and because it allows to define 
a "canonical" BCS order parameter. In fact, as already 
discussed, Eq. ([S]) defining the grand-canonical gap, is 
always zero in the canonical ensemble. Thus following 
Ref. [HI we use as a canonical order parameter 



N 

* = ^2 u a v Q 

a=l 



(62) 



Note that 4" is just half of the concurrence (which is a 
local entanglement measure, see as a review [3J|) which 
has been already calculated with the present method^. 



1 1 




FIG. 7: (Color online) Canonical order parameter $ as a 
function of g for different numbers of pairs N p = N/2. 




FIG. 8: (Color online) Zoom of * in the region 0.15 < g < 
0.25 for several N p = N/2. Inset: Scaling ansatz of reference 
[lllfl7] | and its failure for large N p . 



In the large N limit all these correlators must reduce 
to the value in the grand-canonical ensemble, which from 
Eq. © specialized to e a — a is 



1 



2 y/A 2 + (a- N p ) 2 /N* 



(63) 



where we fixed the chemical potential to /i = N/2 = N p 
and we recall that A is given by Eq. l[55|) . Consequently, 
in the same limit, the canonical order parameter is 



2N„ 



*V„ = 



N p =oo 



lim > u„v a 



NA 



N v —*oo 



dx 



4 7-i VA 2 + (x/2) 2 
N 



NA ^ VI + 4A 2 + 1 _ iVA _ 
~2~ ° S Vl + 4A 2 -T ~~ ~2g~ ~ 4. 9 sinh l/2g 



(64) 



It is evident that \? vanishes when the gap A is zero, 
confirming that in the thermodynamic limit it is a good 
order parameter. 

Our results for u a v a are reported in Figs. [5] and [51 In 
the former each plot consists of the various curves at fixed 
g (=0.1, 0.2, 0.4, 0.7) with varying N p . The latter instead 
shows the g dependence at fixed N p . In Fig. [5] also the 
BCS results for any g are reported for comparison. It 
is evident that for all g the results tend to converge to 
the BCS ones, as they must. However this convergence 
is slower as g is smaller, for example for g = 0.1 the max- 
imum at N — 256 is only 90% of the asymptotic result 
and conversely at g = 0.7 the N p = 16 result is already 
99.8%. These finite N p correlations are symmetric with 
respect to (N + l)/2 by construction. However we point 
out that this will not be true for different level corre- 
lations, while in the grand-canonical ensemble they are 
both symmetric. 

In Fig. [7] we report the order parameter *&/N as a 
function of g for several values of N p , and compare it with 
the BCS result. This figure is exactly the same as the one 
for the concurrence obtained by Dunning et al, 16 , with 



the important difference that they considered only N p < 
34 while we pushed the calculation up to N p — 128. We 
could have calculated these correlations for larger N p , but 
the ones considered are already enough to describe the 
crossover from the mesoscopic to the macroscopic regime. 
In fact, Fig. [7] shows that for N p = 128 'F is almost 
indistinguishable from the BCS one Eq. (|64|) . except for 
very small g that are characterized by the scaling (156[) . 

,iFrom the figure, it is evident that for g < g* ~ 0.18 
the BCS limit is reached from above, whereas for g > g* 
it is approached from below. Exactly at g = g* all the 
curves seem to cross in the same point. This is slightly 
different from what was observed before for a small num- 
ber of particle o 11 ' 17 , where to get a similar crossing the 
order parameter was multiplied by N~ v with 7/ ~ 0.94. 
To clarify this point we zoom in on the crossing point 
in Fig. [H It is evident that for N p > 8 all the curves 
approximately cross in g*, but this is not the case for 
smaller sizes. It is then direct to interpret g* as the value 
of g where the leading finite-size correction of order 1/N 
vanishes (in fact these are clearly negative for large g 
and positive for very small ones). The differences for 
smaller size are due to higher order corrections ~ 1/N 2 . 
This fixed point is thus completely analogous to results 
discussed in the previous section for the ground-state en- 
ergy. A very interesting problem would be to calculate g* 
directly from the finite-size form in an analytical manner 
using the l/Np expansion previously discusse d 22 ' 32 . 

The finite-size scaling \fr ~ found in Refs. [TTIIT7l | 
can clearly not be true for large sizes, since 'F is an exten- 
sive quantity. To check for which sizes it stops working, 
in the inset of Fig. [5] we plot ^N v . All the systems with 
sizes N p < 16 cross indeed at the value of Refs. fTilli"7| 
g CT ~ 0.157, but larger systems clearly deviate from this 
fixed point. We can then safely conclude that this scaling 
ansatz is effective only for N p < 16. In Ref. [ll[ a second 
crossing point has been also found for a larger value of 
the pairing constant. According to our analysis also this 
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FIG. 9: (Color online) Off-diagonal correlations {S 1 as a 
function of a/N. 



FIG. 10: (Color online) (S~ +l S+) as function of a/N. 



fixed point is present only for relatively small number of 
pairs. 

For large g the BCS result is the leading term for large 
N and easily gives ^/N = 1/2 - l/(48.g 2 ) + 0(.g" 3 , iV" 1 ) 
whereas for small coupling we have^ 



N 



ln(3 + v/8)) 



+ 0(l/laN p ), for ff <l. (65) 



Note that for large g we have an N p independent re- 
sult while for small g there is a square-root singularity 
in Np. The latter is again a manifestation of the non- 
perturbative nature of superconductivity. Both these an- 
alytical results are perfectly reproduced by our numerics. 



B. Static correlation functions among different 
levels 

Despite several interesting features of the correlation 
functions among the same levels that we have just dis- 
cussed, these are qualitatively very similar to the grand- 
canonical ones. On the other hand, correlation functions 
between different levels (known as off-diagonal ones) are 
a strong signature of the canonical BCS-like pairing cor- 
relations and should be relevant for the interpretation 
of tunneling experiments. In fact, within the grand- 
canonical ensemble (and so for N = oo) these four-point 
correlation functions factorize to the product of two point 
ones (i.e. in this approximation the Cooper pairs are 
free). Oppositely, in the canonical ensemble they are 
non-trivial functions of both the pairs as a consequence 
of quantum fluctuations. Following Ref. [ll|, we concen- 
trate here on the two correlation functions 



Our results for different values of g and N are reported 
in Fig. l9l and [TTl respectively. 

For N — > oo, as a consequence of factorization, we have 



(S a S%) = U a V a U(3V(3 

1 A 



A 



4 ^A 2 + {a- N p ) 2 /N 2 v/A 2 + {(3- N p f/N 



• (67) 



In particular at fixed /3 these correlations are symmetric 
with respect to a = N p . Including some trivial finite 
size effect, in the grand-canonical ensemble one would 
expect a symmetry at (N + l)/2 as for the on-level corre- 
lations. However, as evident from the figures this is not 
the case in the canonical description. Furthermore, the 
smaller g is the more asymmetrical are the correlations. 
Such asymmetries are very pronounced for all the S~ S + 
correlators, as for example showed in Fig. [TU] where we 
report as the other extreme (compared to {S^S^)), the 
correlator (S^ +1 S^). Thus the asymmetries can be used 
to understand the degree of "canonicality" of a system. 
Note in particular the very different scales in Figs. [H] 
and 1 101 off-diagonal correlation functions are much more 
important when one of the levels is close to the Fermi 
point, a fact that is not surprising being true also in the 
grand-canonical ensemble. 

A last property that is not apparent from the plots but 
that is true (even if not evident from the determinant 
representation) is that 



(S a Sp) — (Sp S&) , 



(68) 



that we checked for all the values we calculated. 



(5f5+), and (S?S*) 



The correlation function S^Sp in the grand-canonical 
(66) ensemble also factorizes into the product of two point 
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FIG. 12: Summed correlation function ^ a p—\ \ S a 

as a function of g for various number of pairs. Inset: small g 

behavior compared with the analytic expression. 



FIG. 11: (Color online) {SfS z a ) as function of a/N. 



ones: 



{a - N p )/2N 



(p-N p )/2N 



N p ) 2 /N 2 



y/A 2 + (a- N p ) 2 /N 2 v'A 2 + Jfi ■ 

(69) 

and it is an odd function at a = M (or (3). Again the 
finite N results do not have this symmetry, that is recov- 
ered only in the thermodynamic limit. The smoothing of 
the step-like structure increasing g is a well-known effect 
also in the grand-canonical ensemble. 



C. Off-diagonal order parameter 

Another fundamental quantity is the so called off- 
diagonal long-range order parameter— defined by 



OD 



1 

N~„ 



N 

E 

a,P=X 



(70) 



that as a difference with 'F takes into account the effect 
of non-diagonal correlations, 'Fod is clearly accessible 
from the direct computation of the off-diagonal correla- 
tion functions (as already proposed^), but this is not 
needed. In fact, it can be obtained without the use 
of the determinant representation, using the Hellmann- 
Feynman theorem (an alternative method of calculation 
has been also proposed—). The derivative of the ground- 
state energy with respect to the coupling strength, allows 
us to directly compute the double sum over all levels of 
the static S + S~ correlation function, i.e. 



9oD = K 



N 

E 



S rv So ) — 



1 dEp(g) 
N P dg 



(71) 



Fig. [12] shows this summed correlation for different 
pair numbers. At g = 0, the only contributing terms 
are the one for which (3 — a since no correlations be- 
tween pairs in different levels exist and we trivially have 

N v . As the interaction is turned on 



inter-level correlations build up rapidly until they satu- 
rate for maximally correlated wavefunctions. For every 
N p , this large g limit is clearly given by ^od = N p + 1. 

The small g behavior can be obtained analytically from 
the known result for the energy^ 

1-Eq = .g + 2 5 2 ln2 + 0( 5 3 ,(lni)" 1 ) 

=> $od = l + <?41n2 + 0(<? 2 ). (72) 

This curve is shown in the inset of Fig. [12] and perfectly 
agrees with the numerical results. Again the deviations 
from this behavior start to occur at a value of g given by 
the usual logarithmic scaling of Eq. ([56]) . 

In the thermodynamic limit, as a consequence of the 
factorization, f ofl is trivially related to 'J as 



OD = 



* N p =oo 



(73) 



signaling that one is extensive if and only if the other one 
is. However, as evident from the figure, for fixed finite 
N p this is not true and the two quantities are indepen- 
dent. Furthermore for small g, in the regime that is not 
accessible by the BCS ansatz, they are both linear in g 
and cannot be in a quadratic relation as for large N p . 
Actually, it has been proven^! that \P and "Foe satisfy 
the following relations for any value of g and N p 



1 

N 



*(* - 1) < *OD < 1 



N 
~Ni 



(74) 



that for N p — > oo are trivial bounds, but not for finite N p 
We checked that our calculations satisfy these bounds. 
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The direct knowledge of ^od allows for a last, very im- 
portant consistency check. In fact the value found from 
Hellmann-Feynman theorem must equal, at half-filling, 

the 

sum TJp Z/a,^=i calculated from the deter- 

minant representation. We checked for all N p < 64 that 
this is indeed the case. 

VI. CONCLUSIONS 

We have studied the static correlation functions of the 
reduced BCS model in the canonical ensemble. From the 
theoretical point of view we simplified the results of Ref. 
[HI giving the correlations as sums over only N r determi- 
nants of N r x N r matrices. This allowed us to calculate 
the correlation functions for very large numbers of par- 
ticles, describing the crossover from mesoscopic to the 
thermodynamic limit, and going beyond previous exact 
or approximate studies. In particular with such accurate 
calculations we were able to discuss critically some con- 
jectured scaling forms for the canonical order parameter. 
For example we rule out the idea of any phase transition 
as a function of the (positive) pairing strength and num- 
ber of particles, in agreement with other analyses based 
on thermodynamical quantities ^ 38 ' 39 . We also calculate 
the off-diagonal long-range order parameter by using the 
Hellmann-Feynman theorem. 

A first interesting step to go beyond what has been 
done here would be to find a single determinant repre- 
sentation for the correlation functions. We made several 
attempts in this direction, but so far unsuccessfully. 

We only analyzed the case of N non-degenerate 



equidistant energy levels at half-filling. This is the 
most interesting model from the condensed matter point 
of view. However for the description of pairing in 
nuclei other choices of the parameters e a are more 
natura l 9 ' 40 ' 41 ' 44 . These could be treated by a simple 
adaptation of our results. 

Furthermore the method presented here allows in prin- 
ciple to obtain dynamical correlation functions as sums 
of form factors over the excited states. This is a more nu- 
merically demanding problem, but can be tackled in the 
same way as other Bethe Ansatz solvable model o 42 ' 43 . 
Also, by summing over the excited states, one can access 
the finite temperature thermodynamics and this can help 
in understanding some open question a 44 ' 45 for ultra-small 
metallic grains. 

The Hamiltonian {T]) is the simplest one with pairing 
terms. More general models with several couplings have 
been proposed to explain complicated pairing in con- 
densed matter— and nuclear physics 4 ^. Some have been 
also shown to be integrable 14,48 . As a consequence it 
would be extremely interesting to tackle these models 
with methods similar to those presented here. 
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